home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / src / X11 / mpeg_play-2.1 / jrevdct.c < prev    next >
C/C++ Source or Header  |  1995-05-09  |  47KB  |  1,622 lines
  1. /*
  2.  * jrevdct.c
  3.  *
  4.  * Copyright (C) 1991, 1992, Thomas G. Lane.
  5.  * This file is part of the Independent JPEG Group's software.
  6.  * For conditions of distribution and use, see the accompanying README file.
  7.  *
  8.  * This file contains the basic inverse-DCT transformation subroutine.
  9.  *
  10.  * This implementation is based on an algorithm described in
  11.  *   C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
  12.  *   Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
  13.  *   Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
  14.  * The primary algorithm described there uses 11 multiplies and 29 adds.
  15.  * We use their alternate method with 12 multiplies and 32 adds.
  16.  * The advantage of this method is that no data path contains more than one
  17.  * multiplication; this allows a very simple and accurate implementation in
  18.  * scaled fixed-point arithmetic, with a minimal number of shifts.
  19.  * 
  20.  * I've made lots of modifications to attempt to take advantage of the
  21.  * sparse nature of the DCT matrices we're getting.  Although the logic
  22.  * is cumbersome, it's straightforward and the resulting code is much
  23.  * faster.
  24.  *
  25.  * A better way to do this would be to pass in the DCT block as a sparse
  26.  * matrix, perhaps with the difference cases encoded.
  27.  */
  28.  
  29. #include <string.h>
  30. #include "video.h"
  31. #include "proto.h"
  32.  
  33. #define GLOBAL            /* a function referenced thru EXTERNs */
  34.   
  35. /* We assume that right shift corresponds to signed division by 2 with
  36.  * rounding towards minus infinity.  This is correct for typical "arithmetic
  37.  * shift" instructions that shift in copies of the sign bit.  But some
  38.  * C compilers implement >> with an unsigned shift.  For these machines you
  39.  * must define RIGHT_SHIFT_IS_UNSIGNED.
  40.  * RIGHT_SHIFT provides a proper signed right shift of an INT32 quantity.
  41.  * It is only applied with constant shift counts.  SHIFT_TEMPS must be
  42.  * included in the variables of any routine using RIGHT_SHIFT.
  43.  */
  44.   
  45. #ifdef RIGHT_SHIFT_IS_UNSIGNED
  46. #define SHIFT_TEMPS    INT32 shift_temp;
  47. #define RIGHT_SHIFT(x,shft)  \
  48.     ((shift_temp = (x)) < 0 ? \
  49.      (shift_temp >> (shft)) | ((~((INT32) 0)) << (32-(shft))) : \
  50.      (shift_temp >> (shft)))
  51. #else
  52. #define SHIFT_TEMPS
  53. #define RIGHT_SHIFT(x,shft)    ((x) >> (shft))
  54. #endif
  55.  
  56. /*
  57.  * This routine is specialized to the case DCTSIZE = 8.
  58.  */
  59.  
  60. #if DCTSIZE != 8
  61.   Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
  62. #endif
  63.  
  64.  
  65. /*
  66.  * A 2-D IDCT can be done by 1-D IDCT on each row followed by 1-D IDCT
  67.  * on each column.  Direct algorithms are also available, but they are
  68.  * much more complex and seem not to be any faster when reduced to code.
  69.  *
  70.  * The poop on this scaling stuff is as follows:
  71.  *
  72.  * Each 1-D IDCT step produces outputs which are a factor of sqrt(N)
  73.  * larger than the true IDCT outputs.  The final outputs are therefore
  74.  * a factor of N larger than desired; since N=8 this can be cured by
  75.  * a simple right shift at the end of the algorithm.  The advantage of
  76.  * this arrangement is that we save two multiplications per 1-D IDCT,
  77.  * because the y0 and y4 inputs need not be divided by sqrt(N).
  78.  *
  79.  * We have to do addition and subtraction of the integer inputs, which
  80.  * is no problem, and multiplication by fractional constants, which is
  81.  * a problem to do in integer arithmetic.  We multiply all the constants
  82.  * by CONST_SCALE and convert them to integer constants (thus retaining
  83.  * CONST_BITS bits of precision in the constants).  After doing a
  84.  * multiplication we have to divide the product by CONST_SCALE, with proper
  85.  * rounding, to produce the correct output.  This division can be done
  86.  * cheaply as a right shift of CONST_BITS bits.  We postpone shifting
  87.  * as long as possible so that partial sums can be added together with
  88.  * full fractional precision.
  89.  *
  90.  * The outputs of the first pass are scaled up by PASS1_BITS bits so that
  91.  * they are represented to better-than-integral precision.  These outputs
  92.  * require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
  93.  * with the recommended scaling.  (To scale up 12-bit sample data further, an
  94.  * intermediate INT32 array would be needed.)
  95.  *
  96.  * To avoid overflow of the 32-bit intermediate results in pass 2, we must
  97.  * have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26.  Error analysis
  98.  * shows that the values given below are the most effective.
  99.  */
  100.  
  101. #ifdef EIGHT_BIT_SAMPLES
  102. #define PASS1_BITS  2
  103. #else
  104. #define PASS1_BITS  1        /* lose a little precision to avoid overflow */
  105. #endif
  106.  
  107. #define ONE    ((INT32) 1)
  108.  
  109. #define CONST_SCALE (ONE << CONST_BITS)
  110.  
  111. /* Convert a positive real constant to an integer scaled by CONST_SCALE.
  112.  * IMPORTANT: if your compiler doesn't do this arithmetic at compile time,
  113.  * you will pay a significant penalty in run time.  In that case, figure
  114.  * the correct integer constant values and insert them by hand.
  115.  */
  116.  
  117. #define FIX(x)    ((INT32) ((x) * CONST_SCALE + 0.5))
  118.  
  119. /* When adding two opposite-signed fixes, the 0.5 cancels */
  120. #define FIX2(x)    ((INT32) ((x) * CONST_SCALE))
  121.  
  122. /* Descale and correctly round an INT32 value that's scaled by N bits.
  123.  * We assume RIGHT_SHIFT rounds towards minus infinity, so adding
  124.  * the fudge factor is correct for either sign of X.
  125.  */
  126.  
  127. #define DESCALE(x,n)  RIGHT_SHIFT((x) + (ONE << ((n)-1)), n)
  128.  
  129. /* Multiply an INT32 variable by an INT32 constant to yield an INT32 result.
  130.  * For 8-bit samples with the recommended scaling, all the variable
  131.  * and constant values involved are no more than 16 bits wide, so a
  132.  * 16x16->32 bit multiply can be used instead of a full 32x32 multiply;
  133.  * this provides a useful speedup on many machines.
  134.  * There is no way to specify a 16x16->32 multiply in portable C, but
  135.  * some C compilers will do the right thing if you provide the correct
  136.  * combination of casts.
  137.  * NB: for 12-bit samples, a full 32-bit multiplication will be needed.
  138.  */
  139.  
  140. #ifdef EIGHT_BIT_SAMPLES
  141. #ifdef SHORTxSHORT_32        /* may work if 'int' is 32 bits */
  142. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT16) (const)))
  143. #endif
  144. #ifdef SHORTxLCONST_32        /* known to work with Microsoft C 6.0 */
  145. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT32) (const)))
  146. #endif
  147. #endif
  148.  
  149. #ifndef MULTIPLY        /* default definition */
  150. #define MULTIPLY(var,const)  ((var) * (const))
  151. #endif
  152.  
  153. #ifndef NO_SPARSE_DCT
  154. #define SPARSE_SCALE_FACTOR  8 
  155. #endif
  156.  
  157. /* Precomputed idct value arrays. */
  158.  
  159. static DCTELEM PreIDCT[64][64];
  160.  
  161.  
  162. /*
  163.  *--------------------------------------------------------------
  164.  *
  165.  * init_pre_idct --
  166.  *
  167.  *  Pre-computes singleton coefficient IDCT values.
  168.  *
  169.  * Results:
  170.  *    None.
  171.  *
  172.  * Side effects:
  173.  *    None.
  174.  *
  175.  *--------------------------------------------------------------
  176.  */
  177. void
  178. init_pre_idct() {
  179.   int i;
  180.  
  181.   for (i=0; i<64; i++) {
  182.     memset((char *) PreIDCT[i], 0, 64*sizeof(DCTELEM));
  183.     PreIDCT[i][i] = 1 << SPARSE_SCALE_FACTOR;
  184.     j_rev_dct(PreIDCT[i]);
  185.   }
  186. }
  187.  
  188. #ifndef NO_SPARSE_DCT
  189.   
  190.  
  191.  
  192. /*
  193.  *--------------------------------------------------------------
  194.  *
  195.  * j_rev_dct_sparse --
  196.  *
  197.  *  Performs the inverse DCT on one block of coefficients.
  198.  *
  199.  * Results:
  200.  *    None.
  201.  *
  202.  * Side effects:
  203.  *    None.
  204.  *
  205.  *--------------------------------------------------------------
  206.  */
  207.  
  208. void
  209. j_rev_dct_sparse (data, pos)
  210.      DCTBLOCK data;
  211.      int pos;
  212. {
  213.   short int val;
  214.   register int *dp;
  215.   register int v;
  216.   int quant;
  217.  
  218. #ifdef SPARSE_AC
  219.   register DCTELEM *dataptr;
  220.   DCTELEM *ndataptr;
  221.   int coeff, rr;
  222.   DCTBLOCK tmpdata, tmp2data;
  223.   DCTELEM *tmpdataptr, *tmp2dataptr;
  224.   int printFlag = 1;
  225. #endif
  226.  
  227.  
  228.   /* If DC Coefficient. */
  229.   
  230.   if (pos == 0) {
  231.     dp = (int *)data;
  232.     v = *data;
  233.     quant = 8;
  234.  
  235.     /* Compute 32 bit value to assign.  This speeds things up a bit */
  236.     if (v < 0) {
  237.         val = -v;
  238.         val += (quant >> 1);
  239.         val /= quant;
  240.         val = -val;
  241.     }
  242.     else {
  243.         val = (v + (quant >> 1)) / quant;
  244.     }
  245.  
  246.     v = ((val & 0xffff) | (val << 16));
  247.  
  248.     dp[0] = v;      dp[1] = v;      dp[2] = v;      dp[3] = v;
  249.     dp[4] = v;      dp[5] = v;      dp[6] = v;      dp[7] = v;
  250.     dp[8] = v;      dp[9] = v;      dp[10] = v;      dp[11] = v;
  251.     dp[12] = v;      dp[13] = v;      dp[14] = v;      dp[15] = v;
  252.     dp[16] = v;      dp[17] = v;      dp[18] = v;      dp[19] = v;
  253.     dp[20] = v;      dp[21] = v;      dp[22] = v;      dp[23] = v;
  254.     dp[24] = v;      dp[25] = v;      dp[26] = v;      dp[27] = v;
  255.     dp[28] = v;      dp[29] = v;      dp[30] = v;      dp[31] = v;
  256.  
  257.     return;
  258.   }
  259.   
  260.   /* Some other coefficient. */
  261.  
  262. #ifdef SPARSE_AC
  263.   dataptr = (DCTELEM *)data;
  264.   coeff = dataptr[pos];
  265.   ndataptr = PreIDCT[pos];
  266.  
  267.   printf (" \n");
  268.   printf ("COEFFICIENT = %3d, POSITION = %2d\n", coeff, pos);
  269.  
  270.   for (v=0; v<64; v++) {
  271.     memcpy((char *) tmpdata, data, 64*sizeof(DCTELEM));
  272.   }
  273.   tmpdataptr = (DCTELEM *)tmpdata;
  274.  
  275.   for (v=0; v<64; v++) {
  276.     memcpy((char *) tmp2data, data, 64*sizeof(DCTELEM));
  277.   }
  278.   tmp2dataptr = (DCTELEM *)tmp2data;
  279.  
  280. #ifdef DEBUG
  281.   printf ("original DCTBLOCK:\n");
  282.   for (rr=0; rr<8; rr++) {
  283.     for (v=0; v<8; v++) {
  284.       if (dataptr[8*rr+v] != tmpdataptr[8*rr+v])
  285.         fprintf(stderr, "Error in copy\n");
  286.       printf ("%3d ", dataptr[8*rr+v]);
  287.     }
  288.     printf("\n");
  289.   }
  290. #endif
  291.  
  292.   printf("\n");
  293.  
  294.   for (rr=0; rr<4; rr++) {
  295.     dataptr[0] = (ndataptr[0] * coeff) >> SPARSE_SCALE_FACTOR;
  296.     dataptr[1] = (ndataptr[1] * coeff) >> SPARSE_SCALE_FACTOR;
  297.     dataptr[2] = (ndataptr[2] * coeff) >> SPARSE_SCALE_FACTOR;
  298.     dataptr[3] = (ndataptr[3] * coeff) >> SPARSE_SCALE_FACTOR;
  299.     dataptr[4] = (ndataptr[4] * coeff) >> SPARSE_SCALE_FACTOR;
  300.     dataptr[5] = (ndataptr[5] * coeff) >> SPARSE_SCALE_FACTOR;
  301.     dataptr[6] = (ndataptr[6] * coeff) >> SPARSE_SCALE_FACTOR;
  302.     dataptr[7] = (ndataptr[7] * coeff) >> SPARSE_SCALE_FACTOR;
  303.     dataptr[8] = (ndataptr[8] * coeff) >> SPARSE_SCALE_FACTOR;
  304.     dataptr[9] = (ndataptr[9] * coeff) >> SPARSE_SCALE_FACTOR;
  305.     dataptr[10] = (ndataptr[10] * coeff) >> SPARSE_SCALE_FACTOR;
  306.     dataptr[11] = (ndataptr[11] * coeff) >> SPARSE_SCALE_FACTOR;
  307.     dataptr[12] = (ndataptr[12] * coeff) >> SPARSE_SCALE_FACTOR;
  308.     dataptr[13] = (ndataptr[13] * coeff) >> SPARSE_SCALE_FACTOR;
  309.     dataptr[14] = (ndataptr[14] * coeff) >> SPARSE_SCALE_FACTOR;
  310.     dataptr[15] = (ndataptr[15] * coeff) >> SPARSE_SCALE_FACTOR;
  311.     dataptr += 16;
  312.     ndataptr += 16;
  313.   }
  314.  
  315.   dataptr = (DCTELEM *) data;
  316.  
  317. #ifdef DEBUG 
  318.   printf ("sparse IDCT:\n");
  319.   for (rr=0; rr<8; rr++) {
  320.     for (v=0; v<8; v++) {
  321.       printf ("%3d ", dataptr[8*rr+v]);
  322.     }
  323.     printf("\n");
  324.   }
  325.   printf("\n");
  326. #endif /* DEBUG */
  327. #else  /* NO_SPARSE_AC */
  328. #ifdef FLOATDCT
  329.   if (qualityFlag)
  330.     float_idct(data);
  331.   else
  332. #endif /* FLOATDCT */
  333.     j_rev_dct(data);
  334. #endif /* SPARSE_AC */
  335.   return;
  336.  
  337. }
  338.  
  339. #else
  340.  
  341. /*
  342.  *--------------------------------------------------------------
  343.  *
  344.  * j_rev_dct_sparse --
  345.  *
  346.  *  Performs the original inverse DCT on one block of 
  347.  *  coefficients.
  348.  *
  349.  * Results:
  350.  *    None.
  351.  *
  352.  * Side effects:
  353.  *    None.
  354.  *
  355.  *--------------------------------------------------------------
  356.  */
  357. void
  358. j_rev_dct_sparse (data, pos) 
  359.      DCTBLOCK data;
  360.      int pos;
  361. {
  362.   j_rev_dct(data);
  363. }
  364. #endif /* SPARSE_DCT */
  365.  
  366.  
  367. #ifndef FIVE_DCT
  368.  
  369. #ifndef ORIG_DCT
  370.  
  371.  
  372. /*
  373.  *--------------------------------------------------------------
  374.  *
  375.  * j_rev_dct --
  376.  *
  377.  *  The inverse DCT function.
  378.  *
  379.  * Results:
  380.  *    None.
  381.  *
  382.  * Side effects:
  383.  *    None.
  384.  *
  385.  *--------------------------------------------------------------
  386.  */
  387. void
  388. j_rev_dct (data)
  389.      DCTBLOCK data;
  390. {
  391.   INT32 tmp0, tmp1, tmp2, tmp3;
  392.   INT32 tmp10, tmp11, tmp12, tmp13;
  393.   INT32 z1, z2, z3, z4, z5;
  394.   INT32 d0, d1, d2, d3, d4, d5, d6, d7;
  395.   register DCTELEM *dataptr;
  396.   int rowctr;
  397.   SHIFT_TEMPS
  398.    
  399.   /* Pass 1: process rows. */
  400.   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  401.   /* furthermore, we scale the results by 2**PASS1_BITS. */
  402.  
  403.   dataptr = data;
  404.  
  405.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  406.     /* Due to quantization, we will usually find that many of the input
  407.      * coefficients are zero, especially the AC terms.  We can exploit this
  408.      * by short-circuiting the IDCT calculation for any row in which all
  409.      * the AC terms are zero.  In that case each output is equal to the
  410.      * DC coefficient (with scale factor as needed).
  411.      * With typical images and quantization tables, half or more of the
  412.      * row DCT calculations can be simplified this way.
  413.      */
  414.  
  415.     register int *idataptr = (int*)dataptr;
  416.     d0 = dataptr[0];
  417.     d1 = dataptr[1];
  418.     if ((d1 == 0) && (idataptr[1] | idataptr[2] | idataptr[3]) == 0) {
  419.       /* AC terms all zero */
  420.       if (d0) {
  421.       /* Compute a 32 bit value to assign. */
  422.       DCTELEM dcval = (DCTELEM) (d0 << PASS1_BITS);
  423.       register int v = (dcval & 0xffff) | (dcval << 16);
  424.       
  425.       idataptr[0] = v;
  426.       idataptr[1] = v;
  427.       idataptr[2] = v;
  428.       idataptr[3] = v;
  429.       }
  430.       
  431.       dataptr += DCTSIZE;    /* advance pointer to next row */
  432.       continue;
  433.     }
  434.     d2 = dataptr[2];
  435.     d3 = dataptr[3];
  436.     d4 = dataptr[4];
  437.     d5 = dataptr[5];
  438.     d6 = dataptr[6];
  439.     d7 = dataptr[7];
  440.  
  441.     /* Even part: reverse the even part of the forward DCT. */
  442.     /* The rotator is sqrt(2)*c(-6). */
  443.     if (d6) {
  444.     if (d4) {
  445.         if (d2) {
  446.         if (d0) {
  447.             /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
  448.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  449.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  450.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  451.  
  452.             tmp0 = (d0 + d4) << CONST_BITS;
  453.             tmp1 = (d0 - d4) << CONST_BITS;
  454.  
  455.             tmp10 = tmp0 + tmp3;
  456.             tmp13 = tmp0 - tmp3;
  457.             tmp11 = tmp1 + tmp2;
  458.             tmp12 = tmp1 - tmp2;
  459.         } else {
  460.             /* d0 == 0, d2 != 0, d4 != 0, d6 != 0 */
  461.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  462.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  463.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  464.  
  465.             tmp0 = d4 << CONST_BITS;
  466.  
  467.             tmp10 = tmp0 + tmp3;
  468.             tmp13 = tmp0 - tmp3;
  469.             tmp11 = tmp2 - tmp0;
  470.             tmp12 = -(tmp0 + tmp2);
  471.         }
  472.         } else {
  473.         if (d0) {
  474.             /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
  475.             tmp2 = MULTIPLY(d6, - FIX2(1.306562965));
  476.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  477.  
  478.             tmp0 = (d0 + d4) << CONST_BITS;
  479.             tmp1 = (d0 - d4) << CONST_BITS;
  480.  
  481.             tmp10 = tmp0 + tmp3;
  482.             tmp13 = tmp0 - tmp3;
  483.             tmp11 = tmp1 + tmp2;
  484.             tmp12 = tmp1 - tmp2;
  485.         } else {
  486.             /* d0 == 0, d2 == 0, d4 != 0, d6 != 0 */
  487.             tmp2 = MULTIPLY(d6, - FIX2(1.306562965));
  488.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  489.  
  490.             tmp0 = d4 << CONST_BITS;
  491.  
  492.             tmp10 = tmp0 + tmp3;
  493.             tmp13 = tmp0 - tmp3;
  494.             tmp11 = tmp2 - tmp0;
  495.             tmp12 = -(tmp0 + tmp2);
  496.         }
  497.         }
  498.     } else {
  499.         if (d2) {
  500.         if (d0) {
  501.             /* d0 != 0, d2 != 0, d4 == 0, d6 != 0 */
  502.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  503.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  504.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  505.  
  506.             tmp0 = d0 << CONST_BITS;
  507.  
  508.             tmp10 = tmp0 + tmp3;
  509.             tmp13 = tmp0 - tmp3;
  510.             tmp11 = tmp0 + tmp2;
  511.             tmp12 = tmp0 - tmp2;
  512.         } else {
  513.             /* d0 == 0, d2 != 0, d4 == 0, d6 != 0 */
  514.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  515.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  516.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  517.  
  518.             tmp10 = tmp3;
  519.             tmp13 = -tmp3;
  520.             tmp11 = tmp2;
  521.             tmp12 = -tmp2;
  522.         }
  523.         } else {
  524.         if (d0) {
  525.             /* d0 != 0, d2 == 0, d4 == 0, d6 != 0 */
  526.             tmp2 = MULTIPLY(d6, - FIX2(1.306562965));
  527.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  528.  
  529.             tmp0 = d0 << CONST_BITS;
  530.  
  531.             tmp10 = tmp0 + tmp3;
  532.             tmp13 = tmp0 - tmp3;
  533.             tmp11 = tmp0 + tmp2;
  534.             tmp12 = tmp0 - tmp2;
  535.         } else {
  536.             /* d0 == 0, d2 == 0, d4 == 0, d6 != 0 */
  537.             tmp2 = MULTIPLY(d6, - FIX2(1.306562965));
  538.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  539.  
  540.             tmp10 = tmp3;
  541.             tmp13 = -tmp3;
  542.             tmp11 = tmp2;
  543.             tmp12 = -tmp2;
  544.         }
  545.         }
  546.     }
  547.     } else {
  548.     if (d4) {
  549.         if (d2) {
  550.         if (d0) {
  551.             /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
  552.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  553.             tmp3 = MULTIPLY(d2, (FIX(1.306562965) + .5));
  554.  
  555.             tmp0 = (d0 + d4) << CONST_BITS;
  556.             tmp1 = (d0 - d4) << CONST_BITS;
  557.  
  558.             tmp10 = tmp0 + tmp3;
  559.             tmp13 = tmp0 - tmp3;
  560.             tmp11 = tmp1 + tmp2;
  561.             tmp12 = tmp1 - tmp2;
  562.         } else {
  563.             /* d0 == 0, d2 != 0, d4 != 0, d6 == 0 */
  564.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  565.             tmp3 = MULTIPLY(d2, (FIX(1.306562965) + .5));
  566.  
  567.             tmp0 = d4 << CONST_BITS;
  568.  
  569.             tmp10 = tmp0 + tmp3;
  570.             tmp13 = tmp0 - tmp3;
  571.             tmp11 = tmp2 - tmp0;
  572.             tmp12 = -(tmp0 + tmp2);
  573.         }
  574.         } else {
  575.         if (d0) {
  576.             /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
  577.             tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
  578.             tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
  579.         } else {
  580.             /* d0 == 0, d2 == 0, d4 != 0, d6 == 0 */
  581.             tmp10 = tmp13 = d4 << CONST_BITS;
  582.             tmp11 = tmp12 = -tmp10;
  583.         }
  584.         }
  585.     } else {
  586.         if (d2) {
  587.         if (d0) {
  588.             /* d0 != 0, d2 != 0, d4 == 0, d6 == 0 */
  589.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  590.             tmp3 = MULTIPLY(d2, (FIX(1.306562965) + .5));
  591.  
  592.             tmp0 = d0 << CONST_BITS;
  593.  
  594.             tmp10 = tmp0 + tmp3;
  595.             tmp13 = tmp0 - tmp3;
  596.             tmp11 = tmp0 + tmp2;
  597.             tmp12 = tmp0 - tmp2;
  598.         } else {
  599.             /* d0 == 0, d2 != 0, d4 == 0, d6 == 0 */
  600.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  601.             tmp3 = MULTIPLY(d2, (FIX(1.306562965) + .5));
  602.  
  603.             tmp10 = tmp3;
  604.             tmp13 = -tmp3;
  605.             tmp11 = tmp2;
  606.             tmp12 = -tmp2;
  607.         }
  608.         } else {
  609.         if (d0) {
  610.             /* d0 != 0, d2 == 0, d4 == 0, d6 == 0 */
  611.             tmp10 = tmp13 = tmp11 = tmp12 = d0 << CONST_BITS;
  612.         } else {
  613.             /* d0 == 0, d2 == 0, d4 == 0, d6 == 0 */
  614.             tmp10 = tmp13 = tmp11 = tmp12 = 0;
  615.         }
  616.         }
  617.     }
  618.     }
  619.  
  620.  
  621.     /* Odd part per figure 8; the matrix is unitary and hence its
  622.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  623.      */
  624.  
  625.     if (d7) {
  626.     if (d5) {
  627.         if (d3) {
  628.         if (d1) {
  629.             /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
  630.             z1 = d7 + d1;
  631.             z2 = d5 + d3;
  632.             z3 = d7 + d3;
  633.             z4 = d5 + d1;
  634.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  635.             
  636.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  637.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  638.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  639.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  640.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  641.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  642.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  643.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  644.             
  645.             z3 += z5;
  646.             z4 += z5;
  647.             
  648.             tmp0 += z1 + z3;
  649.             tmp1 += z2 + z4;
  650.             tmp2 += z2 + z3;
  651.             tmp3 += z1 + z4;
  652.         } else {
  653.             /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
  654.             z2 = d5 + d3;
  655.             z3 = d7 + d3;
  656.             z5 = MULTIPLY(z3 + d5, FIX(1.175875602));
  657.             
  658.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  659.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  660.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  661.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  662.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  663.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  664.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  665.             
  666.             z3 += z5;
  667.             z4 += z5;
  668.             
  669.             tmp0 += z1 + z3;
  670.             tmp1 += z2 + z4;
  671.             tmp2 += z2 + z3;
  672.             tmp3 = z1 + z4;
  673.         }
  674.         } else {
  675.         if (d1) {
  676.             /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
  677.             z1 = d7 + d1;
  678.             z4 = d5 + d1;
  679.             z5 = MULTIPLY(d7 + z4, FIX(1.175875602));
  680.             
  681.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  682.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  683.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  684.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  685.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  686.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  687.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  688.             
  689.             z3 += z5;
  690.             z4 += z5;
  691.             
  692.             tmp0 += z1 + z3;
  693.             tmp1 += z2 + z4;
  694.             tmp2 = z2 + z3;
  695.             tmp3 += z1 + z4;
  696.         } else {
  697.             /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
  698.             z5 = MULTIPLY(d7 + d5, FIX(1.175875602));
  699.  
  700.             tmp0 = MULTIPLY(d7, - FIX2(0.601344887));
  701.             tmp1 = MULTIPLY(d5, - FIX2(0.509795578));
  702.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  703.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  704.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  705.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  706.             
  707.             z3 += z5;
  708.             z4 += z5;
  709.             
  710.             tmp0 += z3;
  711.             tmp1 += z4;
  712.             tmp2 = z2 + z3;
  713.             tmp3 = z1 + z4;
  714.         }
  715.         }
  716.     } else {
  717.         if (d3) {
  718.         if (d1) {
  719.             /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
  720.             z1 = d7 + d1;
  721.             z3 = d7 + d3;
  722.             z5 = MULTIPLY(z3 + d1, FIX(1.175875602));
  723.             
  724.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  725.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  726.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  727.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  728.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  729.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  730.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  731.             
  732.             z3 += z5;
  733.             z4 += z5;
  734.             
  735.             tmp0 += z1 + z3;
  736.             tmp1 = z2 + z4;
  737.             tmp2 += z2 + z3;
  738.             tmp3 += z1 + z4;
  739.         } else {
  740.             /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
  741.             z3 = d7 + d3;
  742.             z5 = MULTIPLY(z3, FIX(1.175875602));
  743.             
  744.             tmp0 = MULTIPLY(d7, - FIX2(0.601344887));
  745.             tmp2 = MULTIPLY(d3, FIX(0.509795579));
  746.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  747.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  748.             z3 = MULTIPLY(z3, - FIX2(0.785694958));
  749.             
  750.             tmp0 += z3;
  751.             tmp1 = z2 + z5;
  752.             tmp2 += z3;
  753.             tmp3 = z1 + z5;
  754.         }
  755.         } else {
  756.         if (d1) {
  757.             /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
  758.             z1 = d7 + d1;
  759.             z5 = MULTIPLY(z1, FIX(1.175875602));
  760.  
  761.             tmp0 = MULTIPLY(d7, - FIX2(1.662939224));
  762.             tmp3 = MULTIPLY(d1, FIX2(1.111140466));
  763.             z1 = MULTIPLY(z1, FIX2(0.275899379));
  764.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  765.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  766.  
  767.             tmp0 += z1;
  768.             tmp1 = z4 + z5;
  769.             tmp2 = z3 + z5;
  770.             tmp3 += z1;
  771.         } else {
  772.             /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
  773.             tmp0 = MULTIPLY(d7, - FIX2(1.387039845));
  774.             tmp1 = MULTIPLY(d7, FIX(1.175875602));
  775.             tmp2 = MULTIPLY(d7, - FIX2(0.785694958));
  776.             tmp3 = MULTIPLY(d7, FIX2(0.275899379));
  777.         }
  778.         }
  779.     }
  780.     } else {
  781.     if (d5) {
  782.         if (d3) {
  783.         if (d1) {
  784.             /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
  785.             z2 = d5 + d3;
  786.             z4 = d5 + d1;
  787.             z5 = MULTIPLY(d3 + z4, FIX(1.175875602));
  788.             
  789.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  790.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  791.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  792.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  793.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  794.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  795.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  796.             
  797.             z3 += z5;
  798.             z4 += z5;
  799.             
  800.             tmp0 = z1 + z3;
  801.             tmp1 += z2 + z4;
  802.             tmp2 += z2 + z3;
  803.             tmp3 += z1 + z4;
  804.         } else {
  805.             /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
  806.             z2 = d5 + d3;
  807.             z5 = MULTIPLY(z2, FIX(1.175875602));
  808.             
  809.             tmp1 = MULTIPLY(d5, FIX2(1.662939225));
  810.             tmp2 = MULTIPLY(d3, FIX2(1.111140466));
  811.             z2 = MULTIPLY(z2, - FIX2(1.387039845));
  812.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  813.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  814.             
  815.             tmp0 = z3 + z5;
  816.             tmp1 += z2;
  817.             tmp2 += z2;
  818.             tmp3 = z4 + z5;
  819.         }
  820.         } else {
  821.         if (d1) {
  822.             /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
  823.             z4 = d5 + d1;
  824.             z5 = MULTIPLY(z4, FIX(1.175875602));
  825.             
  826.             tmp1 = MULTIPLY(d5, - FIX2(0.509795578));
  827.             tmp3 = MULTIPLY(d1, FIX2(0.601344887));
  828.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  829.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  830.             z4 = MULTIPLY(z4, FIX2(0.785694958));
  831.             
  832.             tmp0 = z1 + z5;
  833.             tmp1 += z4;
  834.             tmp2 = z2 + z5;
  835.             tmp3 += z4;
  836.         } else {
  837.             /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
  838.             tmp0 = MULTIPLY(d5, FIX(1.175875602));
  839.             tmp1 = MULTIPLY(d5, FIX2(0.275899380));
  840.             tmp2 = MULTIPLY(d5, - FIX2(1.387039845));
  841.             tmp3 = MULTIPLY(d5, FIX2(0.785694958));
  842.         }
  843.         }
  844.     } else {
  845.         if (d3) {
  846.         if (d1) {
  847.             /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
  848.             z5 = d3 + d1;
  849.  
  850.             tmp2 = MULTIPLY(d3, - FIX(1.451774981));
  851.             tmp3 = MULTIPLY(d1, (FIX(0.211164243) - 1));
  852.             z1 = MULTIPLY(d1, FIX(1.061594337));
  853.             z2 = MULTIPLY(d3, - FIX(2.172734803));
  854.             z4 = MULTIPLY(z5, FIX(0.785694958));
  855.             z5 = MULTIPLY(z5, FIX(1.175875602));
  856.             
  857.             tmp0 = z1 - z4;
  858.             tmp1 = z2 + z4;
  859.             tmp2 += z5;
  860.             tmp3 += z5;
  861.         } else {
  862.             /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
  863.             tmp0 = MULTIPLY(d3, - FIX2(0.785694958));
  864.             tmp1 = MULTIPLY(d3, - FIX2(1.387039845));
  865.             tmp2 = MULTIPLY(d3, - FIX2(0.275899379));
  866.             tmp3 = MULTIPLY(d3, FIX(1.175875602));
  867.         }
  868.         } else {
  869.         if (d1) {
  870.             /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
  871.             tmp0 = MULTIPLY(d1, FIX2(0.275899379));
  872.             tmp1 = MULTIPLY(d1, FIX2(0.785694958));
  873.             tmp2 = MULTIPLY(d1, FIX(1.175875602));
  874.             tmp3 = MULTIPLY(d1, FIX2(1.387039845));
  875.         } else {
  876.             /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
  877.             tmp0 = tmp1 = tmp2 = tmp3 = 0;
  878.         }
  879.         }
  880.     }
  881.     }
  882.  
  883.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  884.  
  885.     dataptr[0] = (DCTELEM) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
  886.     dataptr[7] = (DCTELEM) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
  887.     dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
  888.     dataptr[6] = (DCTELEM) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
  889.     dataptr[2] = (DCTELEM) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
  890.     dataptr[5] = (DCTELEM) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
  891.     dataptr[3] = (DCTELEM) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
  892.     dataptr[4] = (DCTELEM) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
  893.  
  894.     dataptr += DCTSIZE;        /* advance pointer to next row */
  895.   }
  896.  
  897.   /* Pass 2: process columns. */
  898.   /* Note that we must descale the results by a factor of 8 == 2**3, */
  899.   /* and also undo the PASS1_BITS scaling. */
  900.  
  901.   dataptr = data;
  902.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  903.     /* Columns of zeroes can be exploited in the same way as we did with rows.
  904.      * However, the row calculation has created many nonzero AC terms, so the
  905.      * simplification applies less often (typically 5% to 10% of the time).
  906.      * On machines with very fast multiplication, it's possible that the
  907.      * test takes more time than it's worth.  In that case this section
  908.      * may be commented out.
  909.      */
  910.  
  911.     d0 = dataptr[DCTSIZE*0];
  912.     d1 = dataptr[DCTSIZE*1];
  913.     d2 = dataptr[DCTSIZE*2];
  914.     d3 = dataptr[DCTSIZE*3];
  915.     d4 = dataptr[DCTSIZE*4];
  916.     d5 = dataptr[DCTSIZE*5];
  917.     d6 = dataptr[DCTSIZE*6];
  918.     d7 = dataptr[DCTSIZE*7];
  919.  
  920.     /* Even part: reverse the even part of the forward DCT. */
  921.     /* The rotator is sqrt(2)*c(-6). */
  922.     if (d6) {
  923.     if (d4) {
  924.         if (d2) {
  925.         if (d0) {
  926.             /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
  927.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  928.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  929.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  930.  
  931.             tmp0 = (d0 + d4) << CONST_BITS;
  932.             tmp1 = (d0 - d4) << CONST_BITS;
  933.  
  934.             tmp10 = tmp0 + tmp3;
  935.             tmp13 = tmp0 - tmp3;
  936.             tmp11 = tmp1 + tmp2;
  937.             tmp12 = tmp1 - tmp2;
  938.         } else {
  939.             /* d0 == 0, d2 != 0, d4 != 0, d6 != 0 */
  940.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  941.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  942.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  943.  
  944.             tmp0 = d4 << CONST_BITS;
  945.  
  946.             tmp10 = tmp0 + tmp3;
  947.             tmp13 = tmp0 - tmp3;
  948.             tmp11 = tmp2 - tmp0;
  949.             tmp12 = -(tmp0 + tmp2);
  950.         }
  951.         } else {
  952.         if (d0) {
  953.             /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
  954.             tmp2 = MULTIPLY(d6, - FIX2(1.306562965));
  955.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  956.  
  957.             tmp0 = (d0 + d4) << CONST_BITS;
  958.             tmp1 = (d0 - d4) << CONST_BITS;
  959.  
  960.             tmp10 = tmp0 + tmp3;
  961.             tmp13 = tmp0 - tmp3;
  962.             tmp11 = tmp1 + tmp2;
  963.             tmp12 = tmp1 - tmp2;
  964.         } else {
  965.             /* d0 == 0, d2 == 0, d4 != 0, d6 != 0 */
  966.             tmp2 = MULTIPLY(d6, -FIX2(1.306562965));
  967.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  968.  
  969.             tmp0 = d4 << CONST_BITS;
  970.  
  971.             tmp10 = tmp0 + tmp3;
  972.             tmp13 = tmp0 - tmp3;
  973.             tmp11 = tmp2 - tmp0;
  974.             tmp12 = -(tmp0 + tmp2);
  975.         }
  976.         }
  977.     } else {
  978.         if (d2) {
  979.         if (d0) {
  980.             /* d0 != 0, d2 != 0, d4 == 0, d6 != 0 */
  981.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  982.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  983.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  984.  
  985.             tmp0 = d0 << CONST_BITS;
  986.  
  987.             tmp10 = tmp0 + tmp3;
  988.             tmp13 = tmp0 - tmp3;
  989.             tmp11 = tmp0 + tmp2;
  990.             tmp12 = tmp0 - tmp2;
  991.         } else {
  992.             /* d0 == 0, d2 != 0, d4 == 0, d6 != 0 */
  993.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  994.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  995.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  996.  
  997.             tmp10 = tmp3;
  998.             tmp13 = -tmp3;
  999.             tmp11 = tmp2;
  1000.             tmp12 = -tmp2;
  1001.         }
  1002.         } else {
  1003.         if (d0) {
  1004.             /* d0 != 0, d2 == 0, d4 == 0, d6 != 0 */
  1005.             tmp2 = MULTIPLY(d6, - FIX2(1.306562965));
  1006.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  1007.  
  1008.             tmp0 = d0 << CONST_BITS;
  1009.  
  1010.             tmp10 = tmp0 + tmp3;
  1011.             tmp13 = tmp0 - tmp3;
  1012.             tmp11 = tmp0 + tmp2;
  1013.             tmp12 = tmp0 - tmp2;
  1014.         } else {
  1015.             /* d0 == 0, d2 == 0, d4 == 0, d6 != 0 */
  1016.             tmp2 = MULTIPLY(d6, - FIX2(1.306562965));
  1017.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  1018.  
  1019.             tmp10 = tmp3;
  1020.             tmp13 = -tmp3;
  1021.             tmp11 = tmp2;
  1022.             tmp12 = -tmp2;
  1023.         }
  1024.         }
  1025.     }
  1026.     } else {
  1027.     if (d4) {
  1028.         if (d2) {
  1029.         if (d0) {
  1030.             /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
  1031.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  1032.             tmp3 = MULTIPLY(d2, (FIX(1.306562965) + .5));
  1033.  
  1034.             tmp0 = (d0 + d4) << CONST_BITS;
  1035.             tmp1 = (d0 - d4) << CONST_BITS;
  1036.  
  1037.             tmp10 = tmp0 + tmp3;
  1038.             tmp13 = tmp0 - tmp3;
  1039.             tmp11 = tmp1 + tmp2;
  1040.             tmp12 = tmp1 - tmp2;
  1041.         } else {
  1042.             /* d0 == 0, d2 != 0, d4 != 0, d6 == 0 */
  1043.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  1044.             tmp3 = MULTIPLY(d2, (FIX(1.306562965) + .5));
  1045.  
  1046.             tmp0 = d4 << CONST_BITS;
  1047.  
  1048.             tmp10 = tmp0 + tmp3;
  1049.             tmp13 = tmp0 - tmp3;
  1050.             tmp11 = tmp2 - tmp0;
  1051.             tmp12 = -(tmp0 + tmp2);
  1052.         }
  1053.         } else {
  1054.         if (d0) {
  1055.             /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
  1056.             tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
  1057.             tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
  1058.         } else {
  1059.             /* d0 == 0, d2 == 0, d4 != 0, d6 == 0 */
  1060.             tmp10 = tmp13 = d4 << CONST_BITS;
  1061.             tmp11 = tmp12 = -tmp10;
  1062.         }
  1063.         }
  1064.     } else {
  1065.         if (d2) {
  1066.         if (d0) {
  1067.             /* d0 != 0, d2 != 0, d4 == 0, d6 == 0 */
  1068.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  1069.             tmp3 = MULTIPLY(d2, (FIX(1.306562965) + .5));
  1070.  
  1071.             tmp0 = d0 << CONST_BITS;
  1072.  
  1073.             tmp10 = tmp0 + tmp3;
  1074.             tmp13 = tmp0 - tmp3;
  1075.             tmp11 = tmp0 + tmp2;
  1076.             tmp12 = tmp0 - tmp2;
  1077.         } else {
  1078.             /* d0 == 0, d2 != 0, d4 == 0, d6 == 0 */
  1079.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  1080.             tmp3 = MULTIPLY(d2, (FIX(1.306562965) + .5));
  1081.  
  1082.             tmp10 = tmp3;
  1083.             tmp13 = -tmp3;
  1084.             tmp11 = tmp2;
  1085.             tmp12 = -tmp2;
  1086.         }
  1087.         } else {
  1088.         if (d0) {
  1089.             /* d0 != 0, d2 == 0, d4 == 0, d6 == 0 */
  1090.             tmp10 = tmp13 = tmp11 = tmp12 = d0 << CONST_BITS;
  1091.         } else {
  1092.             /* d0 == 0, d2 == 0, d4 == 0, d6 == 0 */
  1093.             tmp10 = tmp13 = tmp11 = tmp12 = 0;
  1094.         }
  1095.         }
  1096.     }
  1097.     }
  1098.  
  1099.     /* Odd part per figure 8; the matrix is unitary and hence its
  1100.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  1101.      */
  1102.     if (d7) {
  1103.     if (d5) {
  1104.         if (d3) {
  1105.         if (d1) {
  1106.             /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
  1107.             z1 = d7 + d1;
  1108.             z2 = d5 + d3;
  1109.             z3 = d7 + d3;
  1110.             z4 = d5 + d1;
  1111.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  1112.             
  1113.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  1114.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  1115.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  1116.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1117.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  1118.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  1119.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  1120.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  1121.             
  1122.             z3 += z5;
  1123.             z4 += z5;
  1124.             
  1125.             tmp0 += z1 + z3;
  1126.             tmp1 += z2 + z4;
  1127.             tmp2 += z2 + z3;
  1128.             tmp3 += z1 + z4;
  1129.         } else {
  1130.             /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
  1131.             z2 = d5 + d3;
  1132.             z3 = d7 + d3;
  1133.             z5 = MULTIPLY(z3 + d5, FIX(1.175875602));
  1134.             
  1135.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  1136.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  1137.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  1138.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  1139.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  1140.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  1141.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  1142.             
  1143.             z3 += z5;
  1144.             z4 += z5;
  1145.             
  1146.             tmp0 += z1 + z3;
  1147.             tmp1 += z2 + z4;
  1148.             tmp2 += z2 + z3;
  1149.             tmp3 = z1 + z4;
  1150.         }
  1151.         } else {
  1152.         if (d1) {
  1153.             /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
  1154.             z1 = d7 + d1;
  1155.             z4 = d5 + d1;
  1156.             z5 = MULTIPLY(d7 + z4, FIX(1.175875602));
  1157.             
  1158.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  1159.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  1160.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1161.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  1162.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  1163.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  1164.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  1165.             
  1166.             z3 += z5;
  1167.             z4 += z5;
  1168.             
  1169.             tmp0 += z1 + z3;
  1170.             tmp1 += z2 + z4;
  1171.             tmp2 = z2 + z3;
  1172.             tmp3 += z1 + z4;
  1173.         } else {
  1174.             /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
  1175.             z5 = MULTIPLY(d5 + d7, FIX(1.175875602));
  1176.  
  1177.             tmp0 = MULTIPLY(d7, - FIX2(0.601344887)); 
  1178.             tmp1 = MULTIPLY(d5, - FIX2(0.509795578));
  1179.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  1180.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  1181.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  1182.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  1183.             
  1184.             z3 += z5;
  1185.             z4 += z5;
  1186.             
  1187.             tmp0 += z3;
  1188.             tmp1 += z4;
  1189.             tmp2 = z2 + z3;
  1190.             tmp3 = z1 + z4;
  1191.         }
  1192.         }
  1193.     } else {
  1194.         if (d3) {
  1195.         if (d1) {
  1196.             /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
  1197.             z1 = d7 + d1;
  1198.             z3 = d7 + d3;
  1199.             z5 = MULTIPLY(z3 + d1, FIX(1.175875602));
  1200.             
  1201.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  1202.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  1203.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1204.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  1205.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  1206.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  1207.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  1208.             
  1209.             z3 += z5;
  1210.             z4 += z5;
  1211.             
  1212.             tmp0 += z1 + z3;
  1213.             tmp1 = z2 + z4;
  1214.             tmp2 += z2 + z3;
  1215.             tmp3 += z1 + z4;
  1216.         } else {
  1217.             /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
  1218.             z3 = d7 + d3;
  1219.             z5 = MULTIPLY(z3, FIX(1.175875602));
  1220.             
  1221.             tmp0 = MULTIPLY(d7, - FIX2(0.601344887)); 
  1222.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  1223.             tmp2 = MULTIPLY(d3, FIX(0.509795579));
  1224.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  1225.             z3 = MULTIPLY(z3, - FIX2(0.785694958));
  1226.             
  1227.             tmp0 += z3;
  1228.             tmp1 = z2 + z5;
  1229.             tmp2 += z3;
  1230.             tmp3 = z1 + z5;
  1231.         }
  1232.         } else {
  1233.         if (d1) {
  1234.             /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
  1235.             z1 = d7 + d1;
  1236.             z5 = MULTIPLY(z1, FIX(1.175875602));
  1237.  
  1238.             tmp0 = MULTIPLY(d7, - FIX2(1.662939224)); 
  1239.             tmp3 = MULTIPLY(d1, FIX2(1.111140466));
  1240.             z1 = MULTIPLY(z1, FIX2(0.275899379));
  1241.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  1242.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  1243.  
  1244.             tmp0 += z1;
  1245.             tmp1 = z4 + z5;
  1246.             tmp2 = z3 + z5;
  1247.             tmp3 += z1;
  1248.         } else {
  1249.             /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
  1250.             tmp0 = MULTIPLY(d7, - FIX2(1.387039845));
  1251.             tmp1 = MULTIPLY(d7, FIX(1.175875602));
  1252.             tmp2 = MULTIPLY(d7, - FIX2(0.785694958));
  1253.             tmp3 = MULTIPLY(d7, FIX2(0.275899379));
  1254.         }
  1255.         }
  1256.     }
  1257.     } else {
  1258.     if (d5) {
  1259.         if (d3) {
  1260.         if (d1) {
  1261.             /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
  1262.             z2 = d5 + d3;
  1263.             z4 = d5 + d1;
  1264.             z5 = MULTIPLY(d3 + z4, FIX(1.175875602));
  1265.             
  1266.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  1267.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  1268.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1269.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  1270.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  1271.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  1272.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  1273.             
  1274.             z3 += z5;
  1275.             z4 += z5;
  1276.             
  1277.             tmp0 = z1 + z3;
  1278.             tmp1 += z2 + z4;
  1279.             tmp2 += z2 + z3;
  1280.             tmp3 += z1 + z4;
  1281.         } else {
  1282.             /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
  1283.             z2 = d5 + d3;
  1284.             z5 = MULTIPLY(z2, FIX(1.175875602));
  1285.  
  1286.             tmp1 = MULTIPLY(d5, FIX2(1.662939225));
  1287.             tmp2 = MULTIPLY(d3, FIX2(1.111140466));
  1288.             z2 = MULTIPLY(z2, - FIX2(1.387039845));
  1289.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  1290.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  1291.             
  1292.             tmp0 = z3 + z5;
  1293.             tmp1 += z2;
  1294.             tmp2 += z2;
  1295.             tmp3 = z4 + z5;
  1296.         }
  1297.         } else {
  1298.         if (d1) {
  1299.             /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
  1300.             z4 = d5 + d1;
  1301.             z5 = MULTIPLY(z4, FIX(1.175875602));
  1302.             
  1303.             tmp1 = MULTIPLY(d5, - FIX2(0.509795578));
  1304.             tmp3 = MULTIPLY(d1, FIX2(0.601344887));
  1305.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  1306.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  1307.             z4 = MULTIPLY(z4, FIX2(0.785694958));
  1308.             
  1309.             tmp0 = z1 + z5;
  1310.             tmp1 += z4;
  1311.             tmp2 = z2 + z5;
  1312.             tmp3 += z4;
  1313.         } else {
  1314.             /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
  1315.             tmp0 = MULTIPLY(d5, FIX(1.175875602));
  1316.             tmp1 = MULTIPLY(d5, FIX2(0.275899380));
  1317.             tmp2 = MULTIPLY(d5, - FIX2(1.387039845));
  1318.             tmp3 = MULTIPLY(d5, FIX2(0.785694958));
  1319.         }
  1320.         }
  1321.     } else {
  1322.         if (d3) {
  1323.         if (d1) {
  1324.             /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
  1325.             z5 = d3 + d1;
  1326.  
  1327.             tmp2 = MULTIPLY(d3, - FIX(1.451774981));
  1328.             tmp3 = MULTIPLY(d1, (FIX(0.211164243) - 1));
  1329.             z1 = MULTIPLY(d1, FIX(1.061594337));
  1330.             z2 = MULTIPLY(d3, - FIX(2.172734803));
  1331.             z4 = MULTIPLY(z5, FIX(0.785694958));
  1332.             z5 = MULTIPLY(z5, FIX(1.175875602));
  1333.             
  1334.             tmp0 = z1 - z4;
  1335.             tmp1 = z2 + z4;
  1336.             tmp2 += z5;
  1337.             tmp3 += z5;
  1338.         } else {
  1339.             /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
  1340.             tmp0 = MULTIPLY(d3, - FIX2(0.785694958));
  1341.             tmp1 = MULTIPLY(d3, - FIX2(1.387039845));
  1342.             tmp2 = MULTIPLY(d3, - FIX2(0.275899379));
  1343.             tmp3 = MULTIPLY(d3, FIX(1.175875602));
  1344.         }
  1345.         } else {
  1346.         if (d1) {
  1347.             /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
  1348.             tmp0 = MULTIPLY(d1, FIX2(0.275899379));
  1349.             tmp1 = MULTIPLY(d1, FIX2(0.785694958));
  1350.             tmp2 = MULTIPLY(d1, FIX(1.175875602));
  1351.             tmp3 = MULTIPLY(d1, FIX2(1.387039845));
  1352.         } else {
  1353.             /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
  1354.             tmp0 = tmp1 = tmp2 = tmp3 = 0;
  1355.         }
  1356.         }
  1357.     }
  1358.     }
  1359.  
  1360.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1361.  
  1362.     dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp3,
  1363.                        CONST_BITS+PASS1_BITS+3);
  1364.     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp10 - tmp3,
  1365.                        CONST_BITS+PASS1_BITS+3);
  1366.     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp11 + tmp2,
  1367.                        CONST_BITS+PASS1_BITS+3);
  1368.     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(tmp11 - tmp2,
  1369.                        CONST_BITS+PASS1_BITS+3);
  1370.     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp12 + tmp1,
  1371.                        CONST_BITS+PASS1_BITS+3);
  1372.     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12 - tmp1,
  1373.                        CONST_BITS+PASS1_BITS+3);
  1374.     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp13 + tmp0,
  1375.                        CONST_BITS+PASS1_BITS+3);
  1376.     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp13 - tmp0,
  1377.                        CONST_BITS+PASS1_BITS+3);
  1378.     
  1379.     dataptr++;            /* advance pointer to next column */
  1380.   }
  1381. }
  1382.  
  1383. #else
  1384.  
  1385.  
  1386.  
  1387. /*
  1388.  *--------------------------------------------------------------
  1389.  *
  1390.  * j_rev_dct --
  1391.  *
  1392.  *  The original inverse DCT function.
  1393.  *
  1394.  * Results:
  1395.  *    None.
  1396.  *
  1397.  * Side effects:
  1398.  *    None.
  1399.  *
  1400.  *--------------------------------------------------------------
  1401.  */
  1402. void
  1403. j_rev_dct (data)
  1404.   DCTBLOCK data;
  1405. {
  1406.   INT32 tmp0, tmp1, tmp2, tmp3;
  1407.   INT32 tmp10, tmp11, tmp12, tmp13;
  1408.   INT32 z1, z2, z3, z4, z5;
  1409.   register DCTELEM *dataptr;
  1410.   int rowctr;
  1411.   SHIFT_TEMPS
  1412.  
  1413.   /* Pass 1: process rows. */
  1414.   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  1415.   /* furthermore, we scale the results by 2**PASS1_BITS. */
  1416.  
  1417.   dataptr = data;
  1418.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  1419.     /* Due to quantization, we will usually find that many of the input
  1420.      * coefficients are zero, especially the AC terms.  We can exploit this
  1421.      * by short-circuiting the IDCT calculation for any row in which all
  1422.      * the AC terms are zero.  In that case each output is equal to the
  1423.      * DC coefficient (with scale factor as needed).
  1424.      * With typical images and quantization tables, half or more of the
  1425.      * row DCT calculations can be simplified this way.
  1426.      */
  1427.  
  1428.     if ((dataptr[1] | dataptr[2] | dataptr[3] | dataptr[4] |
  1429.      dataptr[5] | dataptr[6] | dataptr[7]) == 0) {
  1430.       /* AC terms all zero */
  1431.       DCTELEM dcval = (DCTELEM) (dataptr[0] << PASS1_BITS);
  1432.       
  1433.       dataptr[0] = dcval;
  1434.       dataptr[1] = dcval;
  1435.       dataptr[2] = dcval;
  1436.       dataptr[3] = dcval;
  1437.       dataptr[4] = dcval;
  1438.       dataptr[5] = dcval;
  1439.       dataptr[6] = dcval;
  1440.       dataptr[7] = dcval;
  1441.       
  1442.       dataptr += DCTSIZE;    /* advance pointer to next row */
  1443.       continue;
  1444.     }
  1445.  
  1446.     /* Even part: reverse the even part of the forward DCT. */
  1447.     /* The rotator is sqrt(2)*c(-6). */
  1448.  
  1449.     z2 = (INT32) dataptr[2];
  1450.     z3 = (INT32) dataptr[6];
  1451.  
  1452.     z1 = MULTIPLY(z2 + z3, FIX(0.541196100));
  1453.     tmp2 = z1 + MULTIPLY(z3, - FIX(1.847759065));
  1454.     tmp3 = z1 + MULTIPLY(z2, FIX(0.765366865));
  1455.  
  1456.     tmp0 = ((INT32) dataptr[0] + (INT32) dataptr[4]) << CONST_BITS;
  1457.     tmp1 = ((INT32) dataptr[0] - (INT32) dataptr[4]) << CONST_BITS;
  1458.  
  1459.     tmp10 = tmp0 + tmp3;
  1460.     tmp13 = tmp0 - tmp3;
  1461.     tmp11 = tmp1 + tmp2;
  1462.     tmp12 = tmp1 - tmp2;
  1463.     
  1464.     /* Odd part per figure 8; the matrix is unitary and hence its
  1465.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  1466.      */
  1467.  
  1468.     tmp0 = (INT32) dataptr[7];
  1469.     tmp1 = (INT32) dataptr[5];
  1470.     tmp2 = (INT32) dataptr[3];
  1471.     tmp3 = (INT32) dataptr[1];
  1472.  
  1473.     z1 = tmp0 + tmp3;
  1474.     z2 = tmp1 + tmp2;
  1475.     z3 = tmp0 + tmp2;
  1476.     z4 = tmp1 + tmp3;
  1477.     z5 = MULTIPLY(z3 + z4, FIX(1.175875602)); /* sqrt(2) * c3 */
  1478.     
  1479.     tmp0 = MULTIPLY(tmp0, FIX(0.298631336)); /* sqrt(2) * (-c1+c3+c5-c7) */
  1480.     tmp1 = MULTIPLY(tmp1, FIX(2.053119869)); /* sqrt(2) * ( c1+c3-c5+c7) */
  1481.     tmp2 = MULTIPLY(tmp2, FIX(3.072711026)); /* sqrt(2) * ( c1+c3+c5-c7) */
  1482.     tmp3 = MULTIPLY(tmp3, FIX(1.501321110)); /* sqrt(2) * ( c1+c3-c5-c7) */
  1483.     z1 = MULTIPLY(z1, - FIX(0.899976223)); /* sqrt(2) * (c7-c3) */
  1484.     z2 = MULTIPLY(z2, - FIX(2.562915447)); /* sqrt(2) * (-c1-c3) */
  1485.     z3 = MULTIPLY(z3, - FIX(1.961570560)); /* sqrt(2) * (-c3-c5) */
  1486.     z4 = MULTIPLY(z4, - FIX(0.390180644)); /* sqrt(2) * (c5-c3) */
  1487.     
  1488.     z3 += z5;
  1489.     z4 += z5;
  1490.     
  1491.     tmp0 += z1 + z3;
  1492.     tmp1 += z2 + z4;
  1493.     tmp2 += z2 + z3;
  1494.     tmp3 += z1 + z4;
  1495.  
  1496.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1497.  
  1498.     dataptr[0] = (DCTELEM) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
  1499.     dataptr[7] = (DCTELEM) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
  1500.     dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
  1501.     dataptr[6] = (DCTELEM) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
  1502.     dataptr[2] = (DCTELEM) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
  1503.     dataptr[5] = (DCTELEM) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
  1504.     dataptr[3] = (DCTELEM) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
  1505.     dataptr[4] = (DCTELEM) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
  1506.  
  1507.     dataptr += DCTSIZE;        /* advance pointer to next row */
  1508.   }
  1509.  
  1510.   /* Pass 2: process columns. */
  1511.   /* Note that we must descale the results by a factor of 8 == 2**3, */
  1512.   /* and also undo the PASS1_BITS scaling. */
  1513.  
  1514.   dataptr = data;
  1515.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  1516.     /* Columns of zeroes can be exploited in the same way as we did with rows.
  1517.      * However, the row calculation has created many nonzero AC terms, so the
  1518.      * simplification applies less often (typically 5% to 10% of the time).
  1519.      * On machines with very fast multiplication, it's possible that the
  1520.      * test takes more time than it's worth.  In that case this section
  1521.      * may be commented out.
  1522.      */
  1523.  
  1524. #ifndef NO_ZERO_COLUMN_TEST
  1525.     if ((dataptr[DCTSIZE*1] | dataptr[DCTSIZE*2] | dataptr[DCTSIZE*3] |
  1526.      dataptr[DCTSIZE*4] | dataptr[DCTSIZE*5] | dataptr[DCTSIZE*6] |
  1527.      dataptr[DCTSIZE*7]) == 0) {
  1528.       /* AC terms all zero */
  1529.       DCTELEM dcval = (DCTELEM) DESCALE((INT32) dataptr[0], PASS1_BITS+3);
  1530.       
  1531.       dataptr[DCTSIZE*0] = dcval;
  1532.       dataptr[DCTSIZE*1] = dcval;
  1533.       dataptr[DCTSIZE*2] = dcval;
  1534.       dataptr[DCTSIZE*3] = dcval;
  1535.       dataptr[DCTSIZE*4] = dcval;
  1536.       dataptr[DCTSIZE*5] = dcval;
  1537.       dataptr[DCTSIZE*6] = dcval;
  1538.       dataptr[DCTSIZE*7] = dcval;
  1539.       
  1540.       dataptr++;        /* advance pointer to next column */
  1541.       continue;
  1542.     }
  1543. #endif
  1544.  
  1545.     /* Even part: reverse the even part of the forward DCT. */
  1546.     /* The rotator is sqrt(2)*c(-6). */
  1547.  
  1548.     z2 = (INT32) dataptr[DCTSIZE*2];
  1549.     z3 = (INT32) dataptr[DCTSIZE*6];
  1550.  
  1551.     z1 = MULTIPLY(z2 + z3, FIX(0.541196100));
  1552.     tmp2 = z1 + MULTIPLY(z3, - FIX(1.847759065));
  1553.     tmp3 = z1 + MULTIPLY(z2, FIX(0.765366865));
  1554.  
  1555.     tmp0 = ((INT32) dataptr[DCTSIZE*0] + (INT32) dataptr[DCTSIZE*4]) << CONST_BITS;
  1556.     tmp1 = ((INT32) dataptr[DCTSIZE*0] - (INT32) dataptr[DCTSIZE*4]) << CONST_BITS;
  1557.  
  1558.     tmp10 = tmp0 + tmp3;
  1559.     tmp13 = tmp0 - tmp3;
  1560.     tmp11 = tmp1 + tmp2;
  1561.     tmp12 = tmp1 - tmp2;
  1562.     
  1563.     /* Odd part per figure 8; the matrix is unitary and hence its
  1564.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  1565.      */
  1566.  
  1567.     tmp0 = (INT32) dataptr[DCTSIZE*7];
  1568.     tmp1 = (INT32) dataptr[DCTSIZE*5];
  1569.     tmp2 = (INT32) dataptr[DCTSIZE*3];
  1570.     tmp3 = (INT32) dataptr[DCTSIZE*1];
  1571.  
  1572.     z1 = tmp0 + tmp3;
  1573.     z2 = tmp1 + tmp2;
  1574.     z3 = tmp0 + tmp2;
  1575.     z4 = tmp1 + tmp3;
  1576.     z5 = MULTIPLY(z3 + z4, FIX(1.175875602)); /* sqrt(2) * c3 */
  1577.     
  1578.     tmp0 = MULTIPLY(tmp0, FIX(0.298631336)); /* sqrt(2) * (-c1+c3+c5-c7) */
  1579.     tmp1 = MULTIPLY(tmp1, FIX(2.053119869)); /* sqrt(2) * ( c1+c3-c5+c7) */
  1580.     tmp2 = MULTIPLY(tmp2, FIX(3.072711026)); /* sqrt(2) * ( c1+c3+c5-c7) */
  1581.     tmp3 = MULTIPLY(tmp3, FIX(1.501321110)); /* sqrt(2) * ( c1+c3-c5-c7) */
  1582.     z1 = MULTIPLY(z1, - FIX(0.899976223)); /* sqrt(2) * (c7-c3) */
  1583.     z2 = MULTIPLY(z2, - FIX(2.562915447)); /* sqrt(2) * (-c1-c3) */
  1584.     z3 = MULTIPLY(z3, - FIX(1.961570560)); /* sqrt(2) * (-c3-c5) */
  1585.     z4 = MULTIPLY(z4, - FIX(0.390180644)); /* sqrt(2) * (c5-c3) */
  1586.     
  1587.     z3 += z5;
  1588.     z4 += z5;
  1589.     
  1590.     tmp0 += z1 + z3;
  1591.     tmp1 += z2 + z4;
  1592.     tmp2 += z2 + z3;
  1593.     tmp3 += z1 + z4;
  1594.  
  1595.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1596.  
  1597.     dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp3,
  1598.                        CONST_BITS+PASS1_BITS+3);
  1599.     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp10 - tmp3,
  1600.                        CONST_BITS+PASS1_BITS+3);
  1601.     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp11 + tmp2,
  1602.                        CONST_BITS+PASS1_BITS+3);
  1603.     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(tmp11 - tmp2,
  1604.                        CONST_BITS+PASS1_BITS+3);
  1605.     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp12 + tmp1,
  1606.                        CONST_BITS+PASS1_BITS+3);
  1607.     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12 - tmp1,
  1608.                        CONST_BITS+PASS1_BITS+3);
  1609.     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp13 + tmp0,
  1610.                        CONST_BITS+PASS1_BITS+3);
  1611.     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp13 - tmp0,
  1612.                        CONST_BITS+PASS1_BITS+3);
  1613.     
  1614.     dataptr++;            /* advance pointer to next column */
  1615.   }
  1616. }
  1617.  
  1618.  
  1619. #endif /* ORIG_DCT */
  1620. #endif /* FIVE_DCT */
  1621.  
  1622.